Plotting Observations with Python and Diana

This document aims to show how to use Diana to create observation plots from Python scripts and from within the ipython notebook.

Installation

To get access to the Python modules that expose Diana's plotting features, you need to install the python-diana package from the local package repository. Either use the package manager to do this, or enter the following in a terminal:

apt-get install python-diana

Some examples can also be found in the python-diana-examples package.

Importing Modules

You can use Diana's functionality from Python by importing the bdiana and plotcommands modules from the metno package. If you are using an old version of ipython (earlier than version 1.0) then it is also useful to import the embed function from the metno.ipython_extensions module.



In [1]:

    
# Import the modules needed to use Diana.
from metno import bdiana, plotcommands
# Import an extension that lets us embed an image into a worksheet.
from metno.ipython_extensions import embed

Initialisation and Setup

The bdiana module contains the BDiana class. This provides a convenient interface to Diana's functions and also wraps several initialisation steps into one object. To start using Diana, we simply create an instance of this class.



In [2]:

    
b = bdiana.BDiana()

By default, Diana doesn't know anything about data files or other resources. The locations of data files, definitions of plot types, palettes, and other pieces of information are stored in one or more setup files. We choose one of these and use the setup method to let Diana know where to get data from.



In [3]:

    
b.setup("/etc/diana/setup/diana.setup-COMMON")

The method will return True if the setup file was loaded correctly.

Reading Observations

With the setup file loaded, we can ask for information about observations known to Diana.



In [4]:

    
types = b.getObsTypes()
for i in range(0, len(types), 5):
    print ", ".join(map(lambda n: unicode(n, "latin1"), types[i:i+5]))

Choosing one of these, we can obtain the list of parameters available for it.



In [5]:

    
params = b.getObsParameters("Synop")
for i in range(0, len(types), 5):
    print ", ".join(map(lambda n: unicode(n, "latin1"), params[i:i+5]))

We can also obtain a list of times for which observations are available.



In [6]:

    
times = b.getObsTimes("synop")
times[-3:]

Plotting Observations

As with other plots, we set up a map and a plot area to use with the observations we wish to show.



In [7]:

    
m = plotcommands.Map(map="Gshhs-Auto", backcolour="white", land="on")
m.setOption("land.colour", "flesh")

a = plotcommands.Area(name="Norge")

The observations to plot are described by the Observations plot command. Here, we use the Synop set of observations to show all parameters associated with the observations, and we ask for a certain density with which to show the observations. If we want to show them all then we ask for "allobs".



In [8]:

    
o = plotcommands.Observations()
o.setOption("data", "Synop")
o.setOption("parameter", params)
o.setOption("density", 1.0)
#o.setOption("devfield", True)

As usual, we set the plot commands and the plot time before plotting an image.



In [9]:

    
time = times[-3]

b.setPlotCommands([a.text(), m.text(), o.text()])
b.setPlotTime(time)
im = b.plotImage(500, 500)
embed(im)

Obtaining Observation Positions

It can be useful to obtain the positions of observations that have been prepared for plotting. To do this, we can call the getObsPositions method with the observation type and plot time.



In [10]:

    
pos = b.getObsPositions("Synop", time)

The object returned contains information about the number of observations and their x and y coordinates. The x and y values are in projection coordinates, so we need to transform these using a suitable projection object if we want useful coordinates.



In [11]:

    
proj4string = pos.obsArea.P().getProjDefinition()

import pyproj
pr = pyproj.Proj(proj4string)

The first five values before transformation:



In [12]:

    
zip(pos.xpos[:5], pos.ypos[:5])

The first five values after transformation:



In [13]:

    
v = []
for x, y in zip(pos.xpos[:5], pos.ypos[:5]):
    v.append(pr(x, y, inverse = True))
v

Using Anti-aliasing

Raster image plots can be made to look nicer if anti-aliasing is enabled for observations. Set the antialiasing option to True if you want this. It will make no difference to PDF or SVG plots.



In [14]:

    
o = plotcommands.Observations()
o.setOption("data", "Synop")
o.setOption("parameter", b.getObsParameters("Synop"))
o.setOption("density", 1.0)
o.setOption("antialiasing", True)

l = plotcommands.Label(data="")

We plot it in the usual way, remembering to obtain a suitable time to use for plotting.



In [15]:

    
times = b.getObsTimes("Synop")

a = plotcommands.Area(name = "Europa")
b.setPlotCommands([a.text(), m.text(), o.text(), l.text()])
b.setPlotTime(time)
im = b.plotImage(500, 500)
embed(im)